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OSCILLATIONS  IN  PISTON-DRIVEN  SHEAR  FLOW  OF  A 
NON-NEWTONIAN  FLUID 


David  S.  Malkus 
John  A.  Nohel 
Bradley  J.  Plohr 

Abstract.  In  recent  e.xperiments  cn  piston-driven  shear  flow  of  a  highly  elastic  and  very 
viscous  non-Newtonian  fluid.  Lim  and  Schowalter  observed  nearly  periodic  oscillations  in  the 
particle  velocity  at  the  channel  wall  for  particular  values  of  the  constant  volumetric  flow  rate. 
Such  periodicity  has  been  characterized  as  a  "stick/slip"  phenomenon  caused  by  the  failure 
of  the  fluid  to  adhere  to  the  wall.  We  suggest  an  alternative  explanation  for  these  oscillations 
using  a  dynamic  mathematical  model  based  on  the  Johnson-Segalman-Oldroyd  constituti\e 
relation,  the  key  feature  of  which  is  a  non-monotonic  relationship  between  steady  shear 
stress  and  strain  rate.  The  resulting  three-dimensional  equations  of  motion  and  stress  are 
reduced  to  one  space  dimension,  consistent  with  experimental  results.  In  the  inertialess 
approximation,  the  equations  governing  the  flow  can  be  \iewed  as  a  continuous  family  of 
quadratic  ordinary  differential  equations  coupled  by  the  non-local  constraint  that  fi.xes  the 
volumetric  flow  rate.  Varying  the  flow  rate,  the  numerical  simulation  of  solutions  to  this 
system  exhibits  transitions  to  and  from  a  regime  with  persistent  oscillations  that  compare 
favorably  with  the  Lim-Schowalter  observations.  When  the  time-asymptotic  behavior  is 
cyclic,  large  shear  strain  rates  are  observed  in  a  thin  but  macroscopic  layer  near  the  wall. 
The  transitions  are  explained  using  spectral  analysis  of  the  linear  (infinite  dimensional) 
operator  resulting  from  linearization  of  the  quadratic  system  about  a  discontinuous  steady 
state  with  a  jump  discontinuity  in  the  stress  components.  The  persistent  oscillations  arise 
as  a  Hopf  bifurcation  to  periodic  orbits  as  the  volumetric  flow  rate  is  increased  beyond  a 
critical  value. 


1.  Introduction 


In  their  recent  experiments  on  piston-driven  flow  at  a  fixed  volumetric  flow  rate  of  a  highly 
eleistic  and  very  viscous  non-Newtonian  fluid,  F.  Lim  and  W.  Schowalter  [4,  5]  observed 
persistent,  nearly  periodic  oscillations  in  the  particle  velocity  near  the  channel  wall.  They 
characterized  the  periodicity  as  a  stick/slip  phenomenon;  oscillations  are  thought  to  reflect 
periodic  detachment  and  reattachment  of  the  fluid,  while  the  regular  temporal  variation  of 
normal  stresses  leads  to  spatially  periodic  distortion  of  the  extrudate  [1]. 
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In  this  paper,  our  focus  is  on  mathematical  modeling  of  the  Lim-Schowalter  experiment, 
in  which  the  volumetric  flow  rate  is  prescribed  by  driving  the  fluid  with  a  piston  moving 
at  fixed  speed.  Experimental  data  suggest  that  far  away  from  the  piston,  the  flow  is  essen¬ 
tially  one-dimensional.  Highly  elastic  and  very  viscous  non-Newtonian  materials  are  mod¬ 
eled  mathematically  as  incompressible  viscoelastic  fluids  with  fading  memory  and,  for  the 
present,  under  isothermal  conditions.  The  basis  for  our  explanation  of  persistent  oscillations 
in  piston-driven  flow  is  the  differential  Johnson-Segalman-Oldroyd  constitutive  law,  which  is 
consistent  with  experimental  shear  data;  it  exhibits  a  non-monotonic  relation  between  steady 
shear  stress  and  strain  rate.  This  constitutive  model  has  already  been  used  successfully  by 
us  in  Refs.  [6,  7,  12]  to  provide  an  alternative  explanation  of  “spurt"’  and  related  phenomena, 
originally  observed  in  e.xperiments  by  G.  Vinogradov  et  ai.  [14]  in  pressure-gradient  driven 
flows.  (^^  hile  the  Johnson-Segalman-Oidroyd  constitutive  model  performs  well  in  repro¬ 
ducing  particular  shear  flows  of  highly  elastic  and  viscous  non-Newtonian  fluids,  we  neither 
claim  nor  expect  that  it  will  perform  equally  well  in  orher  flows  (e.g..  extcnsional  flow).) 

Piston-driven  flow  at  a  fixed  volumetric  flow  rate  is  modeled  as  an  instantaneous,  globally 
well-posed  feedback-control  problem  in  one  space  dimension,  the  control  being  the  prescribed 
volumetric  flow  rate  and  the  feedback  being  the  pressure  gradient.  deeper  analytical  un¬ 
derstanding  of  the  governing  equations  is  possible  because  the  system  contains  a  small  pa¬ 
rameter.  the  ratio  of  Reynolds  number  to  Deborah  number.  In  the  inertialess  approximation, 
obtained  by  setting  this  parameter  to  zero,  the  equations  governing  the  flow  can  be  viewed 
as  a  continuous  family  of  quadratic  ordinary  differential  equations  coupled  by  the  non-local 
constraint  that  fixes  the  volumetric  flow  rate;  all  solutions  of  this  system  are  bounded,  even 
for  large,  discontinuous  initial  data.  Numerical  simulations  demonstrate  that  beyond  a  crit¬ 
ical  flow  rate,  the  time-asymptotic  behavior  is  cyclic.  Large  shear  strain  rates  are  observed 
in  a  thin,  but  macroscopic,  layer  near  the  wall;  such  a  layer  also  appears  in  pressure-driven 
spurt  flow. 

Further  analysis  of  the  inertialess  system  links  the  observed  oscillations  to  the  occurrence 
of  a  Hopf  bifurcation  that  spawns  a  periodic  orbit  beyond  a  critical  flow  rate.  To  understand 
this  behavior,  we  observe  that  the  dynamic  equations  admit  discontinuous  steady  states  when 
the  volumetric  flow  rate  is  fixed  (corresponding  to  the  wall  stress  between  the  st^^ady  local 
shear  stress  maximum  and  minimum),  and  we  linearize  the  system  around  a  discontinuous 
steady  state  with  a  single  jump.  Parametrizing  such  discontinuous  steady  states  by  wall 
stress  and  flow  rate,  there  are  regions  in  parameter  space  in  which  the  eigenvalues  of  the 
discrete  spectrum  of  the  linearized  operator  change  from  having  negative  real  parts  to  having 
positive  real  parts,  and  there  is  a  separating  curve  along  which  the  real  part  of  the  eigenvalue 
vanishes.  This  information  is  the  basis  for  showing  in  Ref.  [10]  that  a  Hopf  bifurcation  to  a 
periodic  orbit  occurs  when  the  flow  rate  is  raised  beyond  a  critical  value.  However,  proving 
that  the  resulting  periodic  orbit  is  stable  (i.e.,  that  the  bifurcation  is  supercritical)  remains 
a  challenging  open  problem. 

The  outline  of  this  paper  is  as  follows.  In  Sec.  2,  we  summarize  the  development  the 
mathematical  model  for  unsteady,  isothermal,  piston-driven  planar  channel  flow  of  a  highly 
elastic  and  viscous  non-Newtonian  fluid,  based  on  balance  laws  and  the  Johnson-Segalman- 
Oldroyd  differential  constitutive  law,  with  a  single  relaxation  time,  in  three  dimensions.  The 
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inertialess  approximation  for  such  a  flow  takes  the  form  of  a  one-parameter  quadratic  sys¬ 
tem  of  first-order  functional  differential  equations  in  which  the  fixed  volumetric  flow  rate  is 
imposed  as  the  control  and  the  pressure  gradient  is  the  feedback.  In  Sec.  3.  a  two-parameter 
family  of  discontinuous  steady  state  flows  is  determined.  The  two  parameters  can  be  chosen 
to  be  the  ^•olumetric  flow  rate  and  the  position  of  the  stress  discontinuity.  In  Sec.  4,  we  use  a 
numerical  algorithm  for  the  inertialess  system  to  simulate  piston-driven  flows.  W'e  consider 
the  sudden  start-up  problem,  in  which  the  volumetric  flow  rate  is  instantaneously  raised  from 
zero  and  then  held  fixed.  Upon  varying  the  fixed  volumetric  flow  rate,  the  numerical  simu¬ 
lation  of  solutions  exhibits  transitions  to  and  from  a  regime  with  persistent  oscillations;  this 
phenomenon  compares  favorably  with  the  Lim-Schowalter  observations.  In  Sec.  5.  spectral 
properties  of  the  inertialess  system  linearized  about  particular  discontinuous  steady  states 
are  determined.  The  infinite  dimensional  linearized  operator  is  bounded,  and  we  determine 
criteria  for  its  essential  and  discrete  spectra  analytically.  We  use  a  numerical  method  to 
compute  the  discrete  spectrum  of  the  linearized  operator  from  the  analytical  criterion.  The 
discrete  spectrum  consists  of  a  pair  of  simple,  complex  conjugate  eigenvalues,  and  the  sign 
of  the  real  parts  of  these  eigenvalues  changes  from  negative  to  positive  along  a  curve  within 
the  parameter  plane  for  the  discontinuous  steady  states.  In  Sec.  6.  we  relate  the  results  on 
the  discrete  spectrum  obtained  in  Sec.  5  to  results  of  numerical  simulations  in  Sec.  4.  and  we 
explain  the  experimentally  observed  persistent  oscillations  as  a  Hopf  bifurcation  occurring 
when  the  control  parameter  is  increased  beyond  a  critical  value.  Finally,  in  Sec.  7.  we  use 
the  results  of  Secs.  5  and  6  to  offer  an  alternative  explanation  of  the  mechanism  underlying 
the  Lim-Schowalter  experiment. 

2.  Pistcn-Driven  Shear  Flow 

To  model  the  experiments  of  Lim  and  Schowalter,  we  consider  the  planar  Poiseuille  flow 
a  channel  of  a  highly  elastic  and  very  viscous  incompressible  non-Newtonian  fluid  under 
isothermal  conditions,  satisfying  the  Johnson-Segalman-Oldroyd  constitutive  relations  for 
the  evolution  of  the  extra  stress  tensor.  We  summarize  the  derivation  of  the  system  of 
partial  differential  equations  that  governs  such  a  one-dimensional  shear  flow,  starting  from 
balance  laws  and  constitutive  equations  in  three  dimensions;  for  more  details,  see  Refs.  [6,  7j. 

The  channel  is  aligned  along  the  y-axis  and  extends  between  x  =  —h/2  and  x  =  h/2. 
The  flow  is  assumed  to  be  symmetric  about  the  centerline  x  =  0  of  the  channel,  and  we 
restrict  attention  to  the  interval  x  €  [— /i/2,0].  Since  the  fluid  undergoes  simple  shearing, 
the  velocity  and  stress  variables  are  independent  of  y.  In  particular,  the  velocity  field  is 
V  =  (0.  f(x.f)).  so  that  the  flow  is  incompressible  and  the  conservation  of  mass  equation  is 
automatically  satisfied. 

The  total  shear  stress  on  the  fluid  is  the  sum  of  three  contributions,  an  isotropic  pressure, 
a  Newtonian  stress,  and  an  extra  stress  tt.  The  conservation  of  momentum  in  the  x-direction 
implies  that  the  pressure  takes  the  form  p  =  po(x,  f)  —  /(f)y,  with  /  being  the  pressure 
gradient.  We  adopt  the  Johnson-Segalman-Oldroyd  differential  constitutive  law  with  a  single 
relaxation  time  [13.  2]  to  determine  the  extra  stress.  In  shear  flow,  the  extra  stress  is 
expressible  in  terms  of  two  variables,  the  shear  stress  cr(x,  t)  :=  and  the  principal  normal 
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stress  difference  Z(x,t)  :=  -5(1  -  -  7r“);  here  a  6  (-1, 1)  is  the  slip  parameter 

of  the  fluid.  The  fluid  variables  v,  cr.  and  Z  are  governed  by  the  y-momentum  equation 
and  constitutive  differential  equations.  We  introduce  two  dimensionless  parameters  that 
characterize  the  flow:  a,  the  ratio  of  Reynolds  number  to  Deborah  number;  and  c,  the  ratio 
of  Newtonian  viscosity  to  shear  viscosity.  After  suitable  scaling,  the  governing  equations 
become 


avt  -  ax  =  st’xx  +  /  , 

at  -  {Z  +  l)i-x  =  -a  , 

Zt  +  avx  =  —Z  . 

on  the  interval  [-1/2,0].  The  boundary  conditions  are 

i;(-l/2.t)  =  0  and  t.’x(0.  t)  =  0. 


(JSO) 


{BC) 


At  t  =  0.  we  impose  the  initial  conditions 


i’(x,0)  =  vq{x)  ,  a(j-.O)  =  ao(j)  ,  and  Z{x,0)  =  Zo{x)  .  {1C) 


For  consistency  of  {IC)  with  {BC),  we  assume  that  r'o(— 1/2)  =  0  and  io(0)  =  0;  to  maintain 
continuity  at  x  =  0.  we  require  that  ao(0)  =  0.  This  condition,  together  with  the  boundary 
condition  fj(0,  t)  =  0  and  the  second  equation  of  system  {JSO),  implies  that  a(0.  t)  =  0  for 
all  t.  Thus  symmetry  about  the  centerline  in  maintained. 

In  the  case  of  highly  elastic  polymer  systems,  a  is  very  small,  7  to  10  orders  of  magnitude 
smaller  than  s,  which  itself  is  of  the  order  of  10"^.  Throughout  this  paper,  therefore,  we 
make  the  inertialess  approximation  q  =  0.  This  approximation  has  been  used  successfully 
in  Refs.  [7,  12]  to  study  pressure-driven  flow  in  which  the  pressure  gradient  /  is  prescribed. 
If  a  =  0,  the  momentum  equation  in  system  {JSO),  combined  with  the  boundary  condi¬ 
tions  {BC)  and  a(0,  t)  =  0,  implies  that 


e 


where 

T{x,t):=-f{t)x  (2.2) 

coincides  with  the  total  shear  stress  a  +  eiv-  Using  Eq.  (2.1),  Vx  can  be  eliminated  from 
system  {JSO). 

In  piston-driven  flow,  the  volumetric  flow  rate  Q  is  prescribed,  Q  being  given  by 

Q{t)  :=  2  f  v{xJ)dx  =  -2  /  xVx{x,t)dx  ;  (2.3) 

J-\/2  J-l/2 

here  an  integration  by  parts  was  carried  out  in  the  second  step.  However,  the  pressure 
gradient  /  is  not  known;  instead,  /  adjusts  to  maintain  the  desired  flow  rate  Q.  We  can 
express  /.  or  equivalently  T,  in  terms  of  the  given  flow  rate  Q  as  follows.  Substituting  for 
Vx  in  Eq.  (2.3)  using  Eq.  (2.1)  leads  to  the  feedback  equation 

T{x,  t)  =  -12sQ(0x  +  24j  f°  ya{y,  t)  dy  .  {FB) 

y-i/z 
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where  we  think  of  the  prescribed  Q  as  the  control  and  /  (or  equivalently  T)  as  the  feedback. 
This  equation  is  written  more  simply  as  T  =  Tq  +  Pa,  where  To{x,t)  :=  -12sQ{t)x  and 
where  the  operator  P,  defined  by 

[Pa]  (x,  t]  ;=  24x  /  ya{y,  t)  dy  (2.4) 

J-l/2 

is  a  projection,  in  the  sense  that  P'^a  ;=  P\Pa]  =  Pa. 

Thus  the  system  governing  inertialess  piston-driven  flow  is  the  quadratic  system  of  func¬ 
tional  differential  equations 

a,  =  {Z  +  l)(^^'^-a  , 

T-c^  _  ^  {QFDE) 

T  =  T^  +  Pa, 

together  with  prescribed  the  initial  conditions  (/C)  for  a{x,0),  Z{x.Q). 

Throughout  the  rest  of  this  paper,  we  will  assume  that  the  prescribed  flow  rate  Q{t)  =  Q 
is  a  positive  constant.  It  is  shown  in  Ref.  [9]  that  the  initial  value  problem  for  {QFDE)  is 
globally  well  posed  for  all  time  and  for  initial  data  a^,  Zq  that  are  either  smooth  or  rough 
(in  the  sense  of  being  only  bounded  and  measurable)  and  of  arbitrary  size;  moreover,  the 
stress  components  a,  Z  remain  bounded.  (We  remark  that  the  global  well-posedness  of  the 
full  governing  system  for  piston-driven  flow  with  a  >  0  is  also  discussed  in  Ref.  [9].)  Beyond 
technical  considerations  that  need  not  concern  us  here,  the  key  to  proving  global  existence 
and  boundedness  of  solutions  and  understanding  the  dynamics  is  that  solutions  of  system 
{QFDE)  satisfy  the  Lyapunov-type  identity 

—  +  {Z  -h  1)'|  =  -2  [cr^  +  {Z  +  \  f  -  l]  .  (L) 

Observe  that  Eq.  (Z)  is  independent  of  a,  e,  and  T.  It  follows  easily  from  Eq.  (Z)  that 
a{-,t)  and  Z{-,  t)  remain  bounded  for  as  long  as  the  solution  a,  Z  of  system  {QFDE)  exists. 
This  fact,  together  with  standard  results,  implies  that  every  locally  defined  solution  can 
be  continued  forward  (uniquely)  for  all  time;  moreover,  the  solution  remains  bounded,  the 
bound  depending  only  on  the  size  of  the  initial  data.  However,  although  solutions  to  system 
{QFDE)  exist  globally  and  remain  bounded,  it  is  difficult  to  determine  their  asymptotic 
behavior  of  solutions  as  t  oc.  Indeed,  the  extensive  literature  on  quadratic  differential  and 
functional  differential  equations  does  include  any  results  on  the  global  qualitative  structure 
of  solutions  of  {QFDE).  The  rest  of  this  paper  is  devoted  to  understanding  some  of  the 
interesting  dynamics  of  the  quadratic  system  {QFDE)  by  a  combination  of  analytical  and 
numerical  methods,  and  this  leads  to  an  explanation  of  the  experimental  observations  of  Lim 
and  Schowalter. 

3.  STE.4DY  Flows. 

Our  goal  is  to  study  the  linearized  stability  and  instability  of  steady  state  solutions  of 
system  {QFDE).  It  will  be  clear  that  steady  states  of  {JSO)  and  of  {QFDE)  coincide. 
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Fig.  1:  The  steady  shear  stress  T  =  u'(ii) 
plotted  vs.  shear  strain  rate  Fi. 
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-h/2  X,  position  0 

Fig.  2:  The  velocity  profile  in  steady  flow 
for  the  path  OACW  shown  in  Fig.  1 . 


We  shall  indicate  a  steady  state  solution  using  an  overbar.  Taking  Q{t)  :=  Q  >  0  to  be  a 
constant  and  setting  at  =  0  and  Zj  =  0  in  (QFDE)  shows  that  system  [QFDE)  admits 
steady  states  a,  Z  satisfying  the  following  requirements: 


a  = 


1  +  V 


7t2 


z  = 


vl 


1  +  V 


7t2 


(3.1) 


r  =  -^+£l!.,  T  =  To  +  Pa,  (3.2) 

where  we  define  Tq{x)  :=  —12sQx. 

The  steady  strain  rate  Vx  is  determined  by  solving  the  steady  stress/strain-rate  relation 


w{vx)  =  T  , 


(3.3) 


where  the  function  w  is  defined  by 

iL’(5)  :=  T-f-j +£s  .  (3.4) 

1  + 

When  £  <  1/8,  the  function  w  is  not  monotone,  as  shown  in  Fig.  1.  This  leads  to  the 
occurrence  of  piecewise  smooth  steady  state  solutions  with  simple  jump  discontinuities  in 
Fi,  implying  a  kink  in  the  steady  velocity  F;  such  a  solution  is  indicated  in  Figs.  1  and  2. 
Formulae  (3.1)  imply  that  the  steady  stresses  a  and  Z  are  piecewise  smooth  functions  having 
jump  discontinuities  whenever  the  strain  rate  Fx  does. 

While  there  are  steady  solutions  with  an  arbitrary  number  of  discontinuities,  we  focus 
on  those  with  a  single  discontinuity.  There  is  a  two-parameter  family  of  such  solutions.  It 
proves  convenient  to  take  as  parameters  the  total  stress  T,  :=  T{x.)  at  the  kink  and  the 
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thickness  6  ;=  x,  +  1/2  of  the  layer  of  high  strain  rate.  The  region  of  allowable  values  for 
these  parameters  is  rectangular:  T.  lies  between  Tm  and  Tm,  the  stresses  at  the  minimum 
and  maximum  of  w\  and  6  lies  in  the  interval  [0, 1/2). 

For  each  choice  of  (T,,  6),  the  strain  rate  u,  is  determined  by  solving 

w(v^(x))  =  T,—  (3.5) 

for  each  x  €  [—1/2, 0);  here  x,.  :=  — 1/2+(5,  and  when  there  are  multiple  solutions  of  Eq.  (3.5), 
Vx(x)  is  chosen  on  the  high  or  low  strain-rate  branch  of  w  according  as  x  <  i,  or  x  >  x.. 
The  function  determines  a  and  Z  through  Eqs.  (3.1).  and  the  flow-rate  Q  for  the  solution 
is  calculated  through  Eq.  (2.3). 

A  more  con\  enient  formula  for  Q  can  be  obtained  using  a  change  of  variables  introduced 
by  Yao  [15]:  in  Eq.  (2.3),  change  variables  from  x  to  s  :=  ri(x)  using  the  relation  ~Jx  =  u'(s) 
to  obtain 


where  in  the  present  application 


(3.7) 


Here,  the  interval  of  integration  is  U  ;=  [0,  s^j  U  [s~,  s,,,.],  where  sf  :=  liraj,__j±  rj(x)  (the  low 
and  high  strain-rate  solutions  of  w{$,)  =  T.)  and  :=  Vx{-\)  (the  high  strain-rate  solution 
of  w{su-)  =  T^),  as  in  Fig.  1.  The  integral  in  Eq.  (3.6)  can  be  evaluated  explicitly  [15]: 


^  t  u  (t)n-'(0  dt  =  -  (I  -  £■)  tan*^  ^  J- 


(1  + 


Ic2c2 
3-  « 


To  make  contact  with  experiment,  however,  it  is  preferable  to  use  parameters  diflFerent 
from  r,  and  6,  namely  the  total  stress  Ty,  T,/{1  —  2^)  at  the  channel  wall  and  the  flow¬ 
rate  Q.  (See  Fig.  4a  below.)  The  lines  T,  =  Tm,  —  Tm,  and  5  =  0,  w-hich  bound 
the  allowable  region  in  the  (T.,  5)-plane.  correspond  respectively  to  the  three  curves  labeled 
“bottom-jumping,"  “top-jumping,”  and  “classical  flow”  in  the  {Ty,,  (5)-plane. 


4.  Numerical  Simul.^tion  of  Piston-Driven  Flows. 


In  order  to  motivate  the  analysis  presented  in  subsequent  sections,  we  first  present  re¬ 
sults  of  numerically  simulating  the  Lim-Schowalter  experiment  using  system  (QFDE).  The 
numerical  algorithm  (which  was  derived  in  Ref.  [9])  is  as  follows.  We  predict  Vx  at  time  level 


n  +  1  using  the  formula 


(4.1) 


for  the  k"*  element  in  the  mesh;  here  k  runs  from  1  to  the  number  of  elements.  .V.  and 
is  the  value  of  the  shear  stress  in  the  element  at  time  level  n.  For  the  given  flow  rate  Q, 
the  stress  J*  is  computed  by  evaluating  Pc  in  s3-stem  {QFDE)  using  the  midpoint  rule: 
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i;*=^^^i^(-12jg  +  24'£<{x,+,-x,)(x,  +  l,*,)/2)  .  (4.2) 

We  then  difference  and  ad^•ance  the  stresses  a  and  Z,  using  the  latest  available  value  of 
(t’i)n+l,  by  the  following  “corrector"  scheme: 

f  ^n  +  l  ~  (1  "b  ^t{Zn  +  l)(^i)n+l  i  ^  q\ 

\  Z„+1  =  (1  -  At)Z„  -  Ata„+i(ri)„+i  . 

The  basic  simulation  that  we  perform  using  the  above  algorithm  idealizes  the  instan¬ 
taneous  start-up  of  a  flow.  In  other  words,  the  test  apparatus  is  filled  with  fluid  at  rest, 
and  then  the  flow  rate  is  suddenl}'  increased  from  zero  to  the  prescribed  value  Q.  Therefore 
the  initial  conditions  are  i’o(.r)  =  <to(x)  =  Zq{i)  =  0.  To  compare  our  results  with  the 
e.xperimental  results,  we  compute  the  pressure  gradient  /  for  each  given  flow  rate  Q.  In  the 
simulations  of  {QFDE),  we  observe  four  distinct  flow  regimes,  corresponding  to  different 
values  of  the  prescribed  flow  rate  Q:  the  classical,  spurt  I,  oscillatory,  and  spurt  II  regimes. 
Typical  plots  of  f{t)  vs.  t  in  the  four  regimes  are  shown  in  Figs.  3a-d. 

In  the  classical  regime,  continuous  steady  solutions  are  achieved.  Steady  solutions  are 
also  achieved  in  the  spurt  I  and  II  regimes;  they  are.  however,  discontinuous  and  have  the 
appearance  of  the  spurt  solutions  achieved  in  pressure-gradient  driven  flows  [6,  7].  .A.s  Q  is 
increased  and  the  transition  to  the  oscillatory  regime  is  approached,  the  time  required  for 
the  flow  to  settle  to  a  steady  state  increases,  and  the  settling  is  accompanied  by  damped 
oscillations  in  f{t)  and  in  the  velocity  near  the  wall.  Finally,  at  a  critical  value  of  Q,  the 
oscillations  appear  to  become  undamped  and  fail  to  settle  down;  this  is  the  oscillatory  regime. 
If  Q  is  increased  further,  a  second  transition  is  observed,  and  steady,  spurt-like  solutions  are 
again  achieved;  this  is  the  spurt  II  regime.  (For  technical  reasons.  Fig.  3d  terminates  before 
the  oscillations  have  died  away  to  graphical  accuracy;  however,  the  simulation  has  been 
carried  out  for  much  longer  time,  and  there  is  no  doubt  that  the  oscillations  die  out.) 

We  wish  to  stress  that  because  we  rely  on  numerical  simulation,  and  because  it  is  difiicult 
to  accurately  distinguish  slightly  damped  systems  from  undamped  systems  numerically,  the 
e.\istence  of  a  true  oscillatory  regime  for  the  full  system  {QFDE)  can  only  be  conjectured 
at  this  point. 

We  studied  transitions  between  these  regimes  in  detail,  as  a  function  of  the  imposed 
volumetric  flow  rate  Q,  for  the  polyisoprene  sample  PI-7  of  Ref.  [14]  (i.e.,  e  =  0.001417). 
As  Q  is  increased  from  zero,  the  transition  from  the  classical  regime  to  the  spurt  I  regime 
occurs  at  approximately  Q  =  0.1.  The  transition  between  the  steady  spurt  I  regime  and 
the  oscillatory  regime  occurs  at  a  critical  Q  of  about  0.3,  whereas  the  Uansition  from  the 
oscillatory  regime  back  to  the  steady  spurt  II  regime  occurs  at  a  critical  Q  of  about  1.4.  We 
empheisize  that  Lim  [4]  also  observed  four  separate  flow  regimes  in  his  experiments,  the  third 
of  which  is  oscillatory  and  the  fourth  of  which  is  a  relatively  steady  regime  at  high  shear 
rate  (inferred  from  Q).  Thus  our  numerical  simulation  of  {QFDE)  conforms  qualitatively 
with  experimental  results. 
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5.  Linearized  Stability  Analysis. 


We  consider  a  solution  of  system  (QFDE)  that  is  a  small  perturbation  of  a  piecew’se 
smooth  steady  state  (a,  Z);  for  simplicity,  we  assume  that  this  steady  state  has  a  single  jump 
discontinuity  at  x,,  at  which  point  r(x,)  =  T,.  We  write  a  =  a  +  a  and  Z  =  Z  +  Z.  The 
linearized  equations  for  a  and  Z  are 


;=  J{a.Z,T) 


{Z+l)/e 

-oje 


where 


Jia.Z.T)  — 


-\-{Z  +  \)l£  {T-a)/£ 
{2a-T)/£  -1 


The  operator  L  is  a  bounded  linear  operator  acting  on  the  Hilbert  space  of  pairs  of 
functions  (a.Z)  G  A'  :=  PSi,  x  PS^,.  where  PSi,  :=  H\—l/2.x.)  x  is  the 

space  of  piecewise  smooth  (i.e..  //*)  functions  on  the  interval  [—1/2.0]  with  a  single  jump 
discontinuity  at  x..  The  linearized  stability  of  the  steady  state  is  determined  by  the  spectral 
properties  of  L  which  we  now  summarize;  the  results  stated  below  are  proved  in  Ref.  [10] . 

Since  the  linear  operator  L  is  bounded  on  A*.  L  is  the  generator  of  an  analytic  semigroup, 
and  the  time-asymptotic  behavior  of  solutions  of  the  linear  system  (5.1)  is  determined  by 
the  discrete  spectrum  of  L.  To  find  it,  we  also  need  to  determine  the  essential  spectrum  of  I. 
The  first  result  states  that  the  essential  spectrum  of  I  coincides  with  the  set  of  eigenvalues 
of  the  Jacobian  matrix  J  for  every  x  €  [-1/2,0]. 

Proposition  5.1:  The  essential  spectrum  of  L  is 

OessiL}  =  IJ  I  A  €  C I A  /s  an  eigenvalue  of  J{o{x),  Z(x).  J(x))|  .  (5.3) 

xe(-i/2.oi 

Furthermore,  if  T,  #  Tm  and  T,  Tm  (i.e.,  the  steady  solution  is  neither  bottom  nor  top 
jumping),  then  there  exists  an  rj  >  0  such  that  ReX  <  —rj  for  all  A  G  (XessiL)- 

The  Jacobian  matrix  J  is  precisely  the  matrix  analyzed  in  Ref.  [7].  .According  to  Prop.  3.2 
of  Ref.  [7],  the  eigenvalues  of  J(a(x),  Z(x),T(x))  are  real  and  negative  for  all  x  >  x.,  whereas 
for  X  <  X,  the  eigenvalues  can  be  non-real,  but  always  have  negative  real  parts. 

Next,  we  characterize  the  discrete  spectrum  of  L,  i.e.,  its  isolated  eigenvalues  of  finite 
multiplicity.  .An  element  A  G  <7disc(T)  is  such  that  there  exists  {d,Z)  ^  (0,0)  satisfying 

A(|)=./(^.2,T)(|)  +  (‘^_;)f  )pa.  (5.4) 

Row  manipulations  and  use  of  the  steady  state  relationships  a  =  (Z  -I-  1)(T'  —  ct)/s  and 
Z  =  -a{J  —  a)l£  show  that  this  condition  is  equivalent  to 


A-t-l-f-(Z-|-l)/£’  Z/(7 

(A  -|-  2)(7  (A  -}-  2)Z  -f-  A  A- 


.)(n-( 
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Since  X  <7ess(L),  one  can  solve  Eq.  (5.5)  for  (a,  Z)  in  terms  of  Pa.  The  solution  is 


where 


E(A,x)  := 


Z{X,x)  := 


a\^(  E(A,x) 
Z  l2(A,i:) 


_  [(A  +  2)Z  +  A  +  l](Z  +  l)/g _ 

[(A  +  2)Z  +  A  +  l][A  +  1  +  [Z  +  !)/£■]  —  (A  +  2)Z 

-{X  +  2)a(Z  +  l)/s _ 

[(A  +  2)Z  +  A  +  l][A  +  1  +  [Z  +  1)/^]  —  (A  +  2)Z 


We  emphasize  that  the  denominator  in  these  expressions  vanishes  for  some  x  €  0]  if  and 

only  if  A  6  a^ssiL)-  The  preceding  calculations  lead  to  the  following  characterization  of  the 
discrete  spectrum; 

Lemma  5.2:  A  complex  number  A  belongs  to  adisc(L)  if  and  only  if  A  ^  aessiL)  and 

24  [  x^E(A.x)  dx  =  1  .  (5.8) 

J-i;2 

Moreover,  if  X  6  atiisc(f-),  then  X  has  multiplicity  1. 

The  integral  in  the  eigenvalue  condition  (5.8)  can  be  evaluated  by  changing  variables 
from  X  to  s  :=  Vx{x)  using  the  relation  -/x  =  w{s)  {cf  Eq.  (3.6)): 

1  -  24  ^  x2E( A,  x)  dx  =  ^  y  w{s)^  [l  -  S(A,  s)]  w'is)  ds  ,  (5.9) 


where  the  interval  of  integration  is  U  :=  [0,s^]  U  [s7,Su;]  and 


E(A,s)  := 


_  A  +  1  -  s^ 

+  !)[!  +  s(A  +  1)]  -  [1  -  £•  -  s(A  +  1)2]  s2  +  es^ 


(5.10) 


We  define  $(A,r,,<5)  to  he  ef  /24  times  the  right-hand  side  of  Eq.  (5.9),  i.e., 

*  ''■  Juil+s^)mX  +  l)[l+e{X  +  l)]-[l-€-eiX  +  iy]s^  +  £s^}  ’  ^ 

The  functions  $  depends  on  the  parameters  T,  and  6  through  s~,  s^,  and  s^,,  which  deter¬ 
mine  U.  As  an  immediate  consequence  of  Lemma  5.2,  we  have  the  following  characterization 
of  the  discrete  spectrum  of  L. 

Proposition  5.3;  The  discrete  spectrum  of  L  is 


adisc(L)  =  {  A  e  C\aUL)  I  $(A,  T.,  6)  =  o}  . 


(5.12) 


While  the  integral  in  the  definition  of  can  be  evaluated  explicitly,  its  exact  form  is 
not  illuminating  for  computing  the  discrete  spectrum  analytically,  and  we  determine  it  nu¬ 
merically  as  follows.  A  numerical  scheme  is  used  to  solve  the  equation  $(A,r,,6)  =  0  for 
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and  instability  for  steady  state  solutions 
with  a  single  jump.  The  dots  connected 
by  dashed  lines  represent  data  from  start¬ 
up  simulations  of  system  {QFDE). 


Fig,  4b:  closeup  view  of  the  portion  of 
Fig.  6a  where  start-up  simulations  of  sys¬ 
tem  {QFDE)  leads  to  instability. 


A,  given  T.  and  6.  This  is  a  system  of  two  real  equations,  Re<J>(o  +  bi,T,,6)  =  0  and 
Im  $(a  +  bi.T,,6)  =0,  for  the  real  and  imaginary  parts  of  A  :=  a  -f  bi.  We  accounted  for  the 
possibility  of  having  multiple  roots  by  plotting  the  two  level  curves  Re  $  =  0  and  Im  4'  =  0 
in  the  (c,  6)-plane  using  a  contour  plotter  and  looking  for  their  intersections.  For  the  range 
of  parameters  T,  and  6  we  explored,  there  is  always  exactly  one  intersection,  and  hence  a 
single  root. 

We  are  interested  in  the  possibility  that  this  element  A  of  the  discrete  spectrum  have 
Re  A  >  0,  so  that  the  corresponding  steady  state  solution  (a,  Z)  is  linearly  unstable.  Re¬ 
garding  A  now^  as  a  function  of  T,  and  6,  instability  would  occur  in  a  certain  region  of  the 
(r.,  6)-plane.  The  region  in  of  instability  is  separated  from  the  region  of  stability  by  the  set 
of  (T.,6)  such  that  ReA  =  0.  Thus  we  are  led  to  find  the  curve  of  parameters  (r,,6)  for 
which  there  exists  a  6  €  R  such  that  ^{bi,T,,6)  =  0. 

The  results  of  solving  this  system  of  equations  numerically  is  depicted  in  Fig.  4a.  (The 
viscosity  ratio  is  taken  to  be  s  =  0.001417,  which  corresponds  [3]  to  the  polyisoprene  sam¬ 
ple  PI-7  of  Ref.  [14].)  Instead  of  plotting  the  ReA  =  0  curve  in  the  (T., 6)-plane^  however, 
we  have  used  as  coordinates  the  wall  stress  T^,  ■=■  T»/{1  —  26)  and  the  flow  rate  Q,  as  given 
in  Eq.  (3.6). 


6.  Rel.ation  of  Analysis  to  Simul.^tion. 

To  compare  the  results  from  numerical  simulation  of  (QFDE)  in  Sec.  4  with  the  linearized 
stability  analysis  of  Sec.  5,  we  observe  that  for  each  Q  >  0.  the  numerical  solution  of  a  start¬ 
up  problem  for  (QFDE)  contains  one  parameter.  6,  that  approaches  a  constant  value  as 


OSCILLATIONS  IN  VISCOELASTIC  SHEAR  FLOW 


13 


time  tends  to  infinity.  Associated  with  this  value  of  6  and  of  Q,  there  is  a  unique  steady 
state  (a(x),Z(x))  (with  a  single  jump  discontinuity)  having  layer  thickness  6  and  flow  rate 
Q:  this  steady  state  might  not  be  stable,  in  which  case  the  solution  does  not  approach  this 
steady  state,  but  rather  seems  to  tend  to  a  limit  cycle  (in  the  a,  Z  phase  plane)  that  in 
some  sense  “encircles"  the  unstable  steady  state.  This  steady  state,  be  it  linearly  stable  or 
unstable,  can  be  computed  by  solving  the  equation 

Q  =  Qsteady{T.,6)  (6.1) 

for  r,,  where  Q steady {T,,  8)  is  the  function  defined  by  the  right  side  of  Yao’s  formula  for  the 
flow  rate,  Eq.  (3.6),  in  which  the  relevant  parameters  have  been  written  in  terms  of  T.  and 
6,  using  relations  of  Sec.  3. 

We  have  simulated  a  sequence  of  start-up  problems,  with  increasing  values  of  Q.  and 
computed  T,  for  each  value  of  8  obtained.  This  gives  a  sequence  of  points  in  the  {T».8)- 
plane.  The  results  for  £  =  0.001417  are  plotted  as  the  discrete  points  in  Figs.  4a  and  4b:  the 
equivalent  coordinates  (Tu,..  Q)  have  been  used,  as  described  at  the  end  of  Sec.  5.  Figure  4b 
shows  that  the  curve  defined  by  this  set  of  points  crosses  transversely  into  and  out  of  the 
unstable  region  as  Q  increases.  The  two  crossing  points  correspond  to  the  transition  points 
between  the  oscillatory  regime  and  the  spurt  I  and  II  regimes,  respectively.  Thus  the  nu¬ 
merical  results  suggest  that  the  region  of  linearly  unstable  steady  states  is  explored  by  the 
dynamical  start-up  experiments.  We  emphasize  that  although  some  steady  state  solutions 
of  {QFDE)  are  linearly  unstable,  all  solutions  of  {QFDE)  are  bounded  for  all  time  (see. 
Sec.  2). 

A  plausible  explanation  of  the  observed  behavior,  as  well  as  of  the  correspondence  between  ' 
the  solutions  of  {QFDE)  and  the  linearized  stability  analysis,  is  that  for  increasing  values 
of  the  flow  rate  Q,  {QFDE)  is  undergoing  Hopf  bifurcation  upon  crossing  into  and  out 
of  a  regime  of  periodic  orbits.  We  prove  such  a  result  in  Ref.  [10]  by  showing  how  the 
infinite-dimensional  flow  problem  generated  by  system  {QFDE)  on  the  Hilbert  space  X 
is  reduced  (using  the  center  manifold  theorem)  to  one  for  a  two-dimensional  vector  field 
to  which  the  classical  Hopf  bifurcation  theorem  can  be  applied  (see  [11],  Theorem  1.13). 
However,  confirming  that  these  periodic  orbits  are  stable  (hence  observable  in  a  physical 
experiment),  requires  showing  that  the  bifurcation  to  periodic  orbits  is  supercritical,  and 
this  remains  an  open  problem  for  the  present. 

It  is  interesting  to  note  that  the  frequency  of  these  periodic  orbits  seems  to  have  physical 
significance:  in  Ref.  [9],  the  frequency  is  observed  to  be  proportional  to  1/s;  according  to 
the  dimensional  analysis  leading  to  system  {QFDE)  and  the  assumptions  made  in  fitting 
model  parameters  to  material  data  [3],  this  observation  translates  into  a  prediction  that 
the  dimensional  frequency  of  oscillation  is  independent  of  molecular  weight  in  a  sequence  of 
experiments  varying  molecular  weight  at  the  same  fixed  Q.  This  seems  to  be  the  case  with 
Lim’s  experimental  data  [4]  to  a  reasonable  degree  of  accuracy. 

7.  Conclusion 

^^e  have  presented  a  mathematical  model  aimed  at  explaining  the  experiments  of  Lim 
and  Schowalter.  Our  model  of  piston-driven  channel  shear  flow  of  a  highly  elastic  and  very 
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viscous  non-Newtonian  fluid  uses  a  constitutive  model  characterized  by  a  non-monotonic 
relationship  between  steady  shear  stress  and  strain  rate.  We  have  described  the  reduction 
of  the  three-dimensional  equations  of  motion  and  stress  to  approximating,  one-dimensional 
systems  that  can  be  studied  by  a  combination  of  numerical  and  analytic  methods.  The 
inertialess  approximation  results  in  a  quadratic  system  of  functional  differential  equations 
in  which  the  prescribed  volumetric  flow  rate  is  imposed  by  a  feedback  control.  Numerical 
solution  of  this  system  exhibits  four  distinct  flow  regimes,  as  do  the  experiments  of  Lim 
and  Schowalter.  The  third  of  these  regimes  exhibits  persistent  oscillations  that  compare 
favorably  to  the  Lim-Schowalter  observations. 

If  the  governing  system  (JSO)  and  its  inertial  appro.ximation  (QFDE)  model  the  ob¬ 
servations  of  Lim  and  Schowalter  correctly,  then  the  details  of  the  flow  are  different  from 
what  these  experimentalists  assumed.  Rather  than  a  stick-slip  flow,  the  flow  that  we  predict 
in  the  oscillatory  regime  has  a  thin  “apparent  slip"  layer  that  exists  for  all  time.  The  layer 
is  unstable,  in  the  sense  that  there  are  large  and  persistent  time  variations  in  the  apparent 
slip  velocity;  these  are  associated  with  persistent  oscillations  in  the  pressure  gradient  that  is 
controlled  by  fixing  the  flow  rate  Q.  For  a  certain  range  of  Q,  we  have  provided  convincing 
evidence  that  these  oscillations  are  a  consequence  of  a  Hopf  bifurcation  to  periodic  orbits  by 
a  combination  of  analytical  and  computational  methods.  Moreover,  the  frequency  of  these 
oscillations  agrees  with  observations  of  Lim.  Our  earlier  analysis  of  the  “spurt  phenomenon" 
in  Refs.  [7,  12]  explains  why  the  control  is  the  source  of  the  instability:  the  flow  would  be¬ 
come  steady  if  the  control  were  removed  by  prescribing  the  pressure  gradient  (although,  of 
course,  the  desired  flow  rate  would  then  almost  surely  not  be  achieved). 
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